A fast and reliable method for the calculation of band structure of solids with hybrid 

functionals 
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A simple approximation within the framework of the hybrid methods for the calculation of the 
electronic structure of solids is presented. By considering only the diagonal elements of the pertur- 
bation operator (Hartree-Fock exchange minus semilocal exchange) calculated in the basis of the 
semilocal orbitals, the computational time is drastically reduced, while keeping very well in most 
studied cases the accuracy of the results obtained with hybrid functionals when applied without any 
approximations. 
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Quantum calculations for solids are commonly done 
with the local density approximation (LDA) or the 
generalized gradient approximation (GGA) of density- 
functional theory (DFT). 1,2 While the properties calcu- 
lated using the total energy (e.g., lattice constant) are 
fairly well described by LDA or GGA, the electronic band 
structure can not be described correctly in many cases. 
For instance, it is well known that the semilocal approx- 
imations (LDA and GGA) lead to unoccupied states in 
semiconductors and insulators which are too low in en- 
ergy, and hence, to too small band gap with respect to 
experiment (see, e.g., Ref. [U). More advanced theories 
can provide more accurate band structures The best 
known are the GW approximation to the self energy £ xc 
within the many-body perturbation theory (Ref. |9|) and 
the hybrid functional o 10 ' 11 which consist of a mixing of 
semilocal and Hartree-Fock (i.e., exact) exchange. 

The GW method represents the state-of-the-art^ and 
leads very often to very accurate results in particu- 
lar if it is applied self-consistentlyi^— and with ver- 
tex corrections. 14 However, due to the huge computa- 
tional effort required by a GW calculation, most of the 
GW results from the literature were obtained non-self- 
consistently (i.e., one-shot Go Wo). Within the Go Wo 
method, the quasiparticle energies c G ^° are obtained as 
solutions of the nonlinear equation 
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where ib s \' and e , are the orbitals and the corre- 
sponding energies obtained from a previous DFT calcu- 
lation with a semilocal (SL) functional [in Eq. ((T|), 
ufj = SE^c /Sp] . The nondiagonal terms of the matrix 
of £ xc — v^c are neglected^ The quasiparticle energies 
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are most of the time close to experiment, but there 



are known cases (e.g., NiO 1 ^) were self-consistency is re- 
ally needed. Beside the heavy cost of a GW calcula- 
tion (Go Wo can still be considered as a very expensive 
method), a drawback is that the convergence of the re- 
sults with respect to the number of unoccupied states can 
be extremely slow as, for example, for ZnOj 17 i 18 

The hybrid functionals are becoming more and more 



popular for solids (see, e.g., Ref. [7). On average, the 
accuracy of the electronic band structure they provide is 
quite similar to Go Wo. In hybrid functionals, a fraction 
a x of semilocal exchange is replaced by the Hartree-Fock 
(HF) exchange: 
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The value of a x which leads to the best agreement with 
experiment depends on (a) the system under study, (b) 
the considered property, and (c) the underlying semilo- 
cal functional E^. Most of the time, the value of a x 
lies in the range 0—0.5. Among the best known hybrid 
functionals, there is the so-called PBE O 19 ' 20 (the func- 
tional considered in the present work) , where the semilo- 
cal functional in Eq. © is the GGA of Perdew, Burke, 
and Ernzerhof 21 (PBE) and the amount a x of Hartree- 
Fock exchange is set to 0.25 (Ref. l22l ). Another way of 
constructing hybrid functionals consists of replacing only 
the short-range part of the semilocal exchange by short- 
range Hartree-Fock as proposed in Ref. [23| for the HSE 
functional (also based on PBE), where the short- and 
long-range parts of exchange are defined by splitting the 
Coulomb operator with the error function. Not consider- 
ing the long-range Hartree-Fock exchange avoids techni- 
cal problems and makes the calculations faster. Note that 
the idea of screening the Hartree-Fock exchange was al- 
ready used by Bylander and Kleinman for their screened- 
exchange LDA functional (sX-LDA) , 6 which can be con- 
sidered as a hybrid functional with 100% (a x = 1) of 
short-range Hartree-Fock exchange (see Ref. [U for re- 
cent sX-LDA calculations). However, even without long- 
range Hartree-Fock, the use of hybrid functionals for 
solids leads to calculations which are one or two orders 
of magnitude more expensive than with semilocal func- 
tionals. The tendency of the PBEO functional is to over- 
estimate small (< 3 eV) band gaps and to underestimate 
large (> 10 eV) band gaps (see, e.g., Ref. [25|), and due 
to the neglect of the long-range Hartree-Fock in HSE, 
the HSE band gaps are smaller than the PBEO band 
gaps.— ~— Note that in Ref. H^it was shown that the per- 
formance of PBEO and HSE could be improved by mak- 
ing the fraction a x of Hartree-Fock exchange dependent 
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on either the average of |Vp| j p in the unit cell (as done 
in Ref. |8j for the modified Becke- Johnson potential 2 ^) or 
the static dielectric constant. 

In this work, a fast way of getting the energies of or- 
bitals from hybrid functionals is proposed. In Ref. l29l . 
the implementation of hybrid functionals (screened and 
unscreened) into the WIEN2k code^S which is based on 
the full-potential linearized augmented plane-wave plus 
local orbitals metho d 31 ' 32 to solve the Kohn-Sham equa- 
tions, was reported. The Hartree-Fock method was im- 
plemented following the method of Massidda, Posternak, 
and Baldereschij 33 which is based on the pseudocharge 
method to solve the Poisson equation^ 

The calculation of the matrix of the nonlocal Hartrce- 
Fock operator v* 1 F is done in a second variational proce- 
dure, i.e., the operator a x (v^ F — uf L ) is considered as a 
perturbation and the semilocal orbitals are used as basis 
functions: 

(^|a x (^ F -^ L )|^ k ). (3) 

The second variational procedure, which was also 
adopted for the implementation of the Hartree-Fock 
equations in other LAPW code s 33 ' 35 ' 36 leads to cheaper 
calculations, since in practice the number of orbitals ifjf^ 
which are used for the construction of Eq. ([3]) can be cho- 
sen to be much smaller than the number of LAPW basis 
functions. Since the shape of the orbitals obtained from 
the semilocal and hybrid functionals can be expected to 
be rather similar, the most important matrix elements 
of Eq. © are the diagonal ones. This leads to the sim- 
ple proposition which consists of calculating the orbital 
energies from hybrid functional the following way: 

= £ nk + (^l«x ~ « X SL ) |^>, (4) 

i.e., the nondiagonal terms are neglected, similarly as 
done in Eq. (TTJ) for the G Q W method J£ Obviously, 
within this approximation the orbitals are not updated 
since the matrix of the operator a x (v^ F — fx) [Eq. ©I 
is diagonal, which means that the ip^s are already the 
eigenvectors. In the following, the results obtained with 
Eq. g]) will be named PBE0 (i.e., one-shot PBEO) 
in analogy with GoWq. Compared to a self-consistent 
calculation, using Eq. Q leads to calculations of the 
band structure which are at least two orders of magni- 
tude faster (neglect of the nondiagonal terms and only 
one iteration). Naturally, the PBE orbitals were used as 
the semilocal orbitals in Eq. (QJ. Actually, Eq. (0} is 
also very closely related to the way the exchange part of 
the derivative discontinuity is calculated: 

(5) 

where tp n k and f/Vk' are the orbitals at the valence band 
maximum and conduction band minimum, respectively, 
and w EXX is the multiplicative exact exchange (EXX) 
potential obtained with the optimized effective potential 
method [see Ref. [H for more discussion on Eq. fl5J]. 



TABLE I: Transition energies (in eV) obtained with the PBE, 
PBEO, and PBEOo methods. The experimental value for 
Cu 2 is from Ref. SJ. See Table I of Ref. |H for the other 
solids. 
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In order to test the accuracy of Eq. for the cal- 
culation of the orbital energies, we have considered the 
following semiconductors and insulators (the structure 
and cubic lattice constant are given in parenthesis): Ar 
(fee, 5.260 A), C (diamond, 3.567 A), Si (diamond, 5.430 
A), GaAs (zinc blende, 5.648 A), MgO (rocksalt, 4.207 
A), NaCl (rocksalt, 5.595 A), Cu 2 {PnSm, 4.27 A), 
MnO (Fm3m, 4.445 A), and NiO (Fm3m, 4.171 A). 
The unit cell of CU2O contains six atoms. Formally, 
Cu has a valency of +1 and therefore the Cu-3d shell 
is full, which means that the correlation effects in the 
Cu-3d shell should not play an important role as it is the 
case for CuO. 37 MnO and NiO have a cubic symmetry, 
but by taking into account the antiferromagnetic phase 
(along the [111] direction of the cubic cell), the symmetry 
is reduced to a rhombohedral one (four atoms in the unit 
cell). MnO and NiO are two of the most studied Mott 
insulators, for which the semilocal functionals are very 
inaccurate due to the strongly localized character of the 
3d electrons ^ Detailed descriptions of the structures of 
Cu 2 and MnO/NiO can be found in Refs. H and|U 
respectively. 

The calculated transition energies, along with exper- 
imental values, are shown in Table Q] As expected, the 
PBE transition energies are by far too small compared to 
the experiment, while the PBEO values are much larger 
and closer to the experiment. However, as already men- 
tioned above, PBEO overestimates small transition en- 
ergies (e.g., GaAs) and underestimates large transition 
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FIG. 1: (Color online) Density of states of GaAs calculated 
with the PBE, PBEO, and PBEOo methods. The Fermi energy 
is set at zero. 



FIG. 3: (Color online) Density of states of one spin component 
of NiO calculated with the PBE, PBEO, and PBEOo methods. 
The Fermi energy is set at zero. 
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FIG. 2: (Color online) Density of states of one spin compo- 
nent of MnO calculated with the PBE, PBEO, and PBEOo 
methods. The Fermi energy is set at zero. 



energies (e.g., Ar). Actually, the important observation 
of the present work is that the PBEO transition ener- 
gies calculated with the approximation given by Eq. (|4]) 
(PBEOo in Table HJ are very similar to the self-consistent 
PBEO results. Indeed, the largest difference is for MnO 
(~ 0.2 eV), which, anyway, represents less than 10% of 
the difference between PBE (1.47 eV) and PBEO (4.23 
eV). For the r — > T transition in GaAs and CU2O, the 
difference between PBEO and PBEOo is ^ 0.1 eV and less 



than 0.1 eV in all other cases. In the case of MgO and 
NaCl, the accuracy of Eq. Q is impressive. 

Figures [TJ [2j and [3] show the density of states (DOS) 
of GaAs, MnO, and NiO, respectively. In the case of 
GaAs, the PBEO and PBE0 o DOSs are indistinguishable, 
which can be explained by the fact that already the PBE 
and PBEO DOSs are very similar and differ only by the 
shift of the unoccupied bands. For MnO as well, the 
agreement between PBEO and PBEOo is very good, albeit 
small differences can be seen. In the case of NiO (the 
most difficult case considered in this work), more visible 
differences between PBEO and PBEOo can be observed. 
For instance, the PBEOo DOS in the energy range —7 and 
—3 eV looks intermediate between the PBE and PBEO 
DOSs. There is less O-p states in the latter case. Also, 
the unoccupied Ni-d peaks are higher in energy by ~ 0.8 
eV with PBEO than with PBE0 . 

In Refs. [l6| and [T2I band gaps of 1.1 and 4.8 eV for 
NiO calculated with non-self-consistent [Eq. (TT])] and self- 
consistent GW calculations, respectively, were reported. 
The importance of self-consistency in this case is rather 
extreme. Such a huge difference is not observed for NiO 
with hybrid functionals when using Eq. ((4]) instead of 
doing a self-consistent PBEO calculation. The explana- 
tion lies in the fact that in the case of Eq. only 
the semilocal orbitals are used for the construction of 
the nonlocal Hartree-Fock operator v^ F , while for Go Wo 
[Eq. (jXJ)] both the orbitals and their energies are used 
for the calculation of the self-energy E xc . In the case of 
CU2O it was shown in Ref. [l3|, that with respect to non- 
self-consistent Go W), updating only the orbital energies 
in E xc already increases the band by 0.46 eV (from 1.34 
to 1.80 eV), while updating also the orbitals increases 
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further the band gap, but by a much smaller value (0.17 
eV) . Other such examples can be found in Ref . |42|. How- 
ever, the effect of self-consistency in the GW method can 
be reduced if the input orbitals and energies are calcu- 
lated from the hybrid (see Ref. |4J) or LDA + U (see Ref. 
S3) methods. 

In summary, it has been shown that for calcula- 
tions with hybrid functionals, neglecting the nondiagonal 
terms in the construction of the Hartree-Fock Hamilto- 
nian (within the second variational procedure) leads to a 
very good approximation. By doing so, the orbital ener- 
gies [Eq. Q] are very close to the energies obtained from 
a self-consistent calculation done without any approxima- 
tions. Several types of semiconductors and insulators, in- 
cluding the more difficult cases of antiferromagnetic MnO 
and NiO, have been considered and the accuracy of the 



approximation has been shown to be (very) good in most 
cases. For NiO, the results are less impressive than in the 
other cases. The approximation leads to calculations of 
the band structure of solids with hybrid functionals which 
are at least two orders of magnitude faster than the self- 
consistent calculations. This approximation allows the 
calculation of the electronic band structure of solids with 
hybrid functionals on much larger systems. 
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